Thermodynamics and coherence of a trapped dipolar Fermi gas 
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We develop a meanfield treatment of a polarized trapped Fermi gas with dipole-dipole interactions. Our 
approach is based on self-consistent semiclassical Hartree-Fock theory that accounts for direct and exchange 
interactions. We discuss our procedure for numerically implementing the calculation. We study the thermo- 
dynamic and the first and second order correlation properties of the system. We find that the system entropy 
depends on the trap geometry, allowing the system to be cooled as the trap aspect ratio is increased, and that 
exchange interactions cause the correlation functions to be anisotropic in the low temperature regime. We also 
find that many uniform gas thermodynamic predictions, for which direct interaction effects vanish, are qualita- 
tively unreliable for trapped systems, most notably for oblate traps. We develop a simplified Hartree formalism 
that is applicable to anisotropic harmonic traps. 

PACS numbers: 03.75.Ss, 74.20.Rp, 67.30.H-,05.30.Fk 



I. INTRODUCTION 

There has been phenomenal recent progress towards pro- 
ducing a quantum degenerate Fermi gas of polar molecules 
[1-3]. The long-range nature of the dipolar interaction and 
the strong dipole moments of ground state heteronuclear 
molecules has opened the door to an exciting array of new 
physics [4-11]. Of particular interest is the superconducting 
regime of this system, first considered in Ref. [12] and devel- 
oped in [13-15]. 

Here we are concerned with a Hartree-Fock description of 
the trapped normal system in the semiclassical approximation. 
Such a formalism has been the workhorse theory for comput- 
ing thermodynamic properties in ultra-cold gases with contact 
interactions, however it is considerably more difficult to eval- 
uate in systems with long-range interactions. We begin by 
briefly surveying previous work in this area. Initial work for 
the zero temperature case by Goral et al. [16] used a Thomas- 
Fermi ansatz. In that approach an isotropic momentum dis- 
tribution was assumed causing the exchange (Fock) term to 
vanish, nevertheless the retention of the direct (Hartree) inter- 
action was observed to distort the system density distribution 
from the aspect ratio of the trapping potential. A variational 
treatment by Miyakawa et al. [17] extended the Thomas- 
Fermi ansatz to allow for an ellipsoidal momentum distribu- 
tion and showed that exchange effects tended to distort the 
momentum distribution of the system. Beyond the variational 
approximation, Zhang and Yi [18] solved the full semiclassi- 
cal theory for the case of a cylindrical trap by minimising the 
energy function of the system using a Monte-Carlo method. 
Lin et al. [19] have presented a similar formalism but in terms 
of a self-consistent calculation of the phase-space distribution 
(Wigner function). We also note analytical results by Chan et 
al. [20] developed for a Fermi-liquid description of the dipo- 
lar Fermi gas, and work in the zero-temperature regime on dy- 
namic properties such as collective modes and expansion dy- 
namics [21-24], and sound propagation [25]. There has also 
been some studies using time-dependent Hartree-Fock theory 
for this system [26, 27]. 

For the finite temperature, regime Endo et al. [28] have 
developed a variational approach, formally valid in the non- 



degenerate regime, which allowed them to assess the effect of 
temperature on the position and momentum space distortions. 
Zhang and Yi [29] have also presented results of a full semi- 
classical Hartree-Fock calculation, which they have used to 
explore the system deformation and stability. 



In this paper we present the results of a fully self-consistent 
Hartree-Fock treatment of the trapped gas similar to the ap- 
proach employed in Ref. [29]. As was revealed in [18], 
for many system properties (particularly stability) this level 
of theory provides large corrections over the variational pre- 
dictions. Furthermore, in the (finite temperature) degenerate 
regime T < T F , there are currently no valid variational the- 
ories. This regime is urgently in need of quantitative theo- 
retical analysis as experiments with polar molecules are now 
approaching degeneracy. 



The outline of the paper is as follows. After introducing the 
basic Hartree-Fock formalism we compare our results against 
others in the literature for the distortion of the system position 
and momentum distribution. We then calculate a wide range 
of thermodynamic quantities of the system, including the first 
calculations for entropy and heat capacity of the trapped gas. 
Using these results we show that adiabatic mechanical defor- 
mation of the trapping potential can be used to change the 
temperature of the gas, in particular that compressing the trap 
in the polarization direction will cause the gas to reduce in 
temperature. We then discuss how to define and calculate cor- 
relation functions for the system. Our results show that ex- 
change interaction effects cause the correlation functions to 
be anisotropic in the low temperature regime. We introduce 
a simplified Hartree theory in which exchange effects are ne- 
glected, and present results demonstrating the accuracy and 
applicability of this theory. We give a full account of our nu- 
merical method and the techniques we employ to make the 
calculations tractable and accurate in the Appendix. 



II. THEORY 

A. Formalism 

We consider a gas of spin polarized fermions that interact 
by a long-range dipole-dipole interaction of the form 



U M (x) 



C dd l-3cos 2 8 

4tt |x| 3 



(1) 



where Cdd = d 2 /eo, with d the electric dipole moment, and 
with 6 the angle between the position vector and the polariza- 
tion axis, which we take to be the z direction. The atoms are 
confined within a cylindrically symmetric harmonic trap 



E/(x) 



2^21 



u;(x* + y*)+uiz 



(2) 



with aspect ratio A = oj z /oj p . 

In the semi-classical approach the system at temperature T 
is described by the Wigner function 



W(x,k) 



exp([e(x,k)-^//c B T) + l' 



(3) 



where \x is the chemical potential and the Hartree-Fock dis- 
persion relation is 



e(x,k) 



h 2 k 2 
2m 



t/(x) + $ D (x)-$ B (x,k), (4) 



with 



<Mx) = J ^"^(x - x ')^(x, k'), 
#u(x, k) = / — -p U dd (k - k')W(x, k'), 



(5) 
(6) 



the direct and exchange interaction terms, respectively. In 
Eq. (6) Udd ls me Fourier transform of the dipole-dipole in- 
teraction, given by 



tf dd (k) = ^[3cos 2 k -l], 



(7) 



where 6^ is the angle between k and k z . 

To find equilibrium solutions Eqs. (3)-(6) must be solved 
self-consistently subject to the additional constraint of atom 
number, i.e. 



N 



dxdk 
— 3^(x,k), 



(2tt) 3 



(8) 



fixed by adjusting the chemical potential. 

We note that the Hartree-Fock dispersion relation (4) can 
be derived by minimizing the free energy (e.g. see [30]) for 
which the total energy is given by [16] 



E = 



dx.dk 
J2nf 



rtfk 2 



2m 



C/(x) 



$u(x) $e(x,k) 



W(x,k). 
(9) 



B. Numerical treatment 

The semiclassical approach has seen extensive application 
to ultra-cold Bose and Fermi gases with contact interactions, 
and has become the de facto standard theory for providing 
thermodynamic information. The numerical solution of the 
semi-classical theory for the dipolar gas is rather more diffi- 
cult, with the following main challenges: 

(i) Interactions are non-local requiring careful choice of 
numerical grids and techniques to deal with the required 
convolutions. 

(ii) The exchange interaction term ($£;) depends on both 
position and momentum variables. 

In regard to point (i), extensive work on dipolar Bose and 
Fermi gases (e.g. see [31-33]) has seen the development of 
numerical tools for dealing with the direct potential (<&£>), 
such as the use of Hankel transforms for exploiting the cylin- 
drical symmetry and allowing convolutions to be performed 
in Fourier space efficiently (e.g. see [34]). However, point (ii) 
remains a considerable challenge and demands a self consis- 
tent solution of W(x, k) in full x- and k-space for the trapped 
case. In contrast, for the case of fermions or bosons with 
local interactions $_e(x, k) —} $e(x) and the momentum 
dependence of the excitations is only through the kinetic en- 
ergy term. Since this momentum dependence is of a simple 
(isotropic) form it can be integrated out to provide a descrip- 
tion only dependent on the density (see Sec. IIIC). Thus, the 
full Wigner function does not need to be constructed to obtain 
a self-consistent solution. Furthermore, with local interactions 
the density is only a function of the local value of the external 
potential and the calculation can be mapped to a single scalar 
variable avoiding the need for a full spatial representation. 

In the dipolar gas, symmetry is broken by the polarization 
direction of the dipoles and in general the most symmetry 
the system can exhibit (even in an isotropic harmonic trap) 
is cylindrical. Because the symmetry axis of the trap we con- 
sider (2) coincides with the polarization direction we can ex- 
ploit this symmetry to simplify the Wigner function to a func- 
tion four variables {p, z, k p , k z }, i.e. cylindrical coordinates 
in x and k space. 

In implementing our numerical algorithm we make con- 
siderable use of Fourier transform techniques and quadrature 
based on discrete cosine and Bessel functions. Extensive de- 
tails of the numerical algorithm are given in Appendix A. 



III. RESULTS 

In this section we present our results for the thermody- 
namic properties of a trapped dipolar Fermi gas. Follow- 
ing the standard convention (e.g. see [29]) we characterize 
our interaction strength in terms of dimensionless parameter 
D t = iV 1 ' 6 Cdd/47r/ki;a 3 10 , where u> is the geometric mean 
trap frequency and a} lo = -Jh/mui. Another energy scale 
we use is the ideal Fermi temperature T F = Tuo^N) 1 / 3 /ks 



for the harmonically trapped gas. By parameterizing the in- 
teraction in terms of D t , the temperature in units of T F and 
lengths in units of A rl / 3 a^, , it is easy to show that the analysis 
is independent of N (see [16]). 

The interplay between trapping geometry and interactions 
is of central importance in the dipolar gas. To explore this re- 
lationship we present results for the cases of prolate (A = 0.1), 
spherical (A = 1) and oblate (A = 10) trapping. As the dipo- 
lar interaction has an attractive component, the system is not 
stable for all parameter regimes, as has been studied exten- 
sively (e.g. see [16-18]). In general, oblate trapping allows 
the largest values of Dt before the onset of collapse, as this 
geometry enhances the repulsive aspect of the long-range in- 
teractions. It has also been shown that thermal fluctuations 
tend to stabilize the system somewhat against collapse [29], 
so that the strictest condition on stability occurs at T = 0. For 
values of D t < 1, all three trapping geometries we consider 
are stable at T = 0. For D t = 2 only the oblate configuration 
of A = 10 is stable at T = 0, while the A = 1 and A = 0.1 



cases become stable at about 0.3Tp 



0.5T F . 



A. Thermal properties of a trapped dipolar Fermi gas 

1. Position and momentum space distortion 

The anisotropy of the dipole-dipole interaction is reflected 
in the equilibrium distribution of the system. This is conve- 
niently characterized by two simple parameters (see [29]) 



(10) 

(11) 



quantifying the momentum and position space distortion of 
the system, respectively, where 




dx.dk 



x 2 W(x,k), 



(12) 



is the x variance, etc. The trap anisotropy appears in the def- 
inition of /3 so that in the absence of interactions, where the 
position density has the anisotropy imposed by the trap, we 
have /3 = 1. In contrast the non-interacting momentum distri- 
bution is spherically symmetric. 

In Fig. 1 we show our results for the deformation param- 
eters for the same parameter regime considered in Fig. 5 of 
Ref. [29]. We find our results are visually identical to their 
results and a quantitative comparison of data reveals a maxi- 
mum relative difference of 0.11%. Our results show that ther- 
mal fluctuations almost completely wash out the interaction 
induced deformations as the Fermi temperature is approached. 
We also find good agreement with the variation predictions for 
deformation by Miyakawa et al. [17] at T = and by Endo et 
al. [28] near T = T F (shown in Fig. 1). As the finite temper- 
ature variational result is based on a Boltzmann description 
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T/T° F 



FIG. 1. (Color online) Deformation of (a) momentum and (b) 
position density distributions. Calculation parameters: A = 0.1 
(blue/dark grey lines), 1 (green/light grey lines) and 10 (red/grey 
lines). The dimensionless interaction strength is D t = 1. Thin 
curves for T > OATp use the variational approach of [28] and 
crosses at T = use the variational approach of [17]. 



it is only valid in the high temperature regime T 3> T F and 
provides a useful check of the asymptotic behavior. We note 
that, at the variational minimum in [28], the kinetic energy is 
always equal to the high temperature limit, 3NksT/2, which 
may affect the accuracy of predictions for a. 



2. Direct and exchange energy 

For many of the results presented in this paper we observe 
an appreciable difference in the behavior of thermodynamic 
parameters of the trapped system from the uniform gas (see 
[29]). The underlying reason is the rather different role of 
interactions in the two systems: For the uniform system the 
direct interaction term vanishes because the spatial density is 



uniform. In contrast, in the trapped system both direct and 
exchange terms contribute. To demonstrate their importance it 
is useful to define their individual contributions to the system 
energy as [see Eq. (9)] 



„ If dxdk , . , 

^ = 2 J /(2^ $D(xWx ' k) ' 



1 f dx.dk 



(2tt): 



$ £ (x,k)^(x,k), 



(13) 
(14) 



which we shall refer to as the direct and exchange energies, 
respectively. 



3 

> 



3 

fe; 
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T/T° F 



by the trap geometry and is always negative. Except for nearly 
spherical traps, the magnitude of Eg tends to be much smaller 
than Ed- 

To explain these observations we begin by considering the 
direct energy. The sign and strength of this quantity is con- 
trolled by the trap geometry, through its influence on the sys- 
tem spatial density profile: A prolate spatial density causes 
the attractive component of Udd (i.e. interactions with dipoles 
in a head-to-tail configuration, see Fig. 3(a)) to dominate and 
Ed is negative, and increasing in magnitude as the system 
becomes more prolate. An oblate spatial density causes the 
repulsive component of Udd (i-e. interactions with dipoles in a 
side-by-side configuration, see Fig. 3(c)) to dominate and Ed 
is positive, and increasing as the system becomes more oblate. 
Note that the spatial density profile has a geometry that is usu- 
ally close to the trap aspect ratio A, however interactions cause 
additional distortion (as has already characterized by the pa- 
rameter /3). For example, in the spherical trap interactions will 
cause the spatial density to deform to be slightly prolate and 
Ed will be slightly negative, as we see in Fig 2. 

A rather similar geometric argument can be applied to the 
exchange interaction, but now for the momentum space dis- 
tribution: It is energetically favorable for the momentum dis- 
tribution to compress along k p and expand along k z . Interest- 
ingly for bosons, where the negative sign does not accompany 
the exchange term, the opposite momentum behavior would 
be expected. However, any anisotropy in the momentum den- 
sity arises solely from the dipole interaction itself in all cases 
we consider 1 and \Ee\ tends to be much smaller than \Ed\, 
except in nearly spherical traps. 



FIG. 2. (Color online) Direct and exchange energy versus tempera- 
ture for a dipolar Fermi gas. Direct energy Ed (solid) and exchange 
energy Ee (dashed). Aspect ratios A = 0.1 (blue/dark grey lines), 1 
(green/light grey lines) and 10 (red/grey lines). Interaction parameter 
Dt — 1. Note that the Ed and Ee results for the A = 1 case are 
almost indistinguishable. 



(a) prolate 



(b) spherical 



(c) oblate 





FIG. 3. (Color online) Schematic showing how the distribution of 
dipoles in (a) prolate, (b) spherical, and (c) oblate geometries. 



3. Chemical potential 

We present results for the low temperature chemical poten- 
tial of the trapped dipolar gas in Table I for the main system 
parameters considered in this paper. We consider the ratio 
n/kBTp, which gives us a quantitative measure of the effect 
of interactions on the Fermi temperature of the interacting sys- 
tem. In practice we evaluate these results at the small but finite 
temperature of T = O.OlTp. We do this because the lowest 
temperature we can solve for is limited by the ability of the 
computational grids we use to resolve the sharp Fermi surface. 
However, as shown in Fig. 4, the chemical rapidly saturates as 
T ->• and the values calculated at T = 0.01T£ should be 
very close to the T = value. 

The behavior of the chemical potential has been considered 
for the uniform gas in [29]. In that work it was found that in- 
creasing the dipolar interaction strength suppressed /j,. In the 
trapped system we observe that prolate trapping geometry can 
enhance this suppression, while an oblate geometry can in- 



We present results for these energies in Fig 2. We ob- 
serve that the direct interaction energy is strongly effected by 
the trap geometry, is significantly increased in magnitude in 
highly anisotropic traps and can be both positive and nega- 
tive. The exchange interaction energy is only slightly affected 



1 The isotropy of mass ensures the non-interacting momentum distribution 
is isotropic, however this could be changed, e.g. by using an optical lattice 
to modify the effective mass in different directions. It may be possible to 
extend our meanfield analysis to this case, e.g. see [35]. 





A 




D t 


0.1 1 


10 


0.5 


0.95 1.00 


1.07 


1 


0.87 0.98 


1.14 


2 


N/A N/A 


1.24 



TABLE I. Low temperature chemical potential fi/ksTp for the sys- 
tem parameters considered in this paper. Chemical potential is eval- 
uated from our Hartree-Fock solution at a temperature of T = 



0.017$. 




0.4 0.6 

T/T° 

FIG. 4. (Color online) Chemical potential (/u/fesTp) as a function of 
temperature. Calculation parameters: A = 0.1 (blue/dark grey lines), 
1 (green/light grey lines) and 10 (red/grey lines). The dimensionless 
interaction strength is D t — 1. Values at T — 0.017$ shown with 



stead cause the chemical potential to increase with increasing 
interaction strength. 



4. Entropy: Mechanism for mechanical cooling 
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FIG. 5. (Color online) Entropy versus temperature for a dipolar 
Fermi gas. (a) Entropy versus temperature, (b) Difference in en- 
tropy from the ideal Fermi gas. Inset shows a magnification of the 
entropy curves with a adiabatic process from initial condition (i) to 
final condition (f) indicated. Aspect ratios A = 0.1 (blue/dark grey 
lines), 1 (green/light grey lines) and 10 (red/grey lines); Interaction 
parameter D f = 0.5 (dotted), D t — 1 (solid), D t — 2 (dashed). 



It is of interest to know how the entropy of the system de- 
pends on the other parameters, such as the trapping potential 
and interaction strength. We can directly calculate the entropy 
from the Wigner function as 

S = -J^{w(x,k)lnW(x,k) 

+ [1-W r (x,k)]ln[l-W(x,k)]}. (15) 

In Fig. 5(a) we show results for the entropy versus tempera- 
ture for three different trap geometries with constant interac- 
tion strength. In Ref. [29] the entropic behavior of the uni- 
form gas was studied and in their numerical simulations they 
found that as the strength of the dipolar interactions were in- 
creased the value of entropy decreased (at fixed temperature). 



We find even richer behavior in the trapped system. When the 
position space density is prolate (i.e. A < 1) interactions tend 
to decrease the entropy (most easily seen by considering the 
difference in entropy form the non-interacting case shown in 
Fig. 5(b)). In contrast, for an oblate density distribution the 
behavior is the opposite of the uniform case and the entropy 
increases with increasing interaction strength. 

We have also computed the cases shown in Fig. 5(b) but 
neglecting exchange interactions (see Sec. IIIC), and find 
that universally this increases the value of entropy by a small 
amount (relative to results including exchange). This demon- 
strates that the entropic shifts we observe are dominated by 
direct interactions, consistent with the energetic observations 
made in Sec. Ill A 2. 



Interestingly, for all cases shown in Fig. 5(b) the differ- 
ence arising from interactions is maximum at a temperature 
of T « 0.2T F , with little dependence on D t or A. This co- 
incides with the temperature at which we observe the most 
rapid change in the deformation parameters (see Fig. 1). Thus, 
0.2T F appears to set the appropriate temperature scale for ex- 
periments to achieve in order to observe the onset of strong 
dipolar effects in thermodynamic properties. Currently this is 
approximately an order of magnitude colder than the coldest 
reported temperature for polar molecules. 

We observe the interesting feature that as the trap ratio A in- 
creases the entropy curves shift to lower temperatures. Thus, 
if the trap geometry is changed adiabatically (constant en- 
tropy) then the system temperature should decrease (see pro- 
cess indicated in inset to Fig. 5(a)). This realizes a mechanical 
cooling process for the system induced by squeezing the dipo- 
lar gas along the polarization direction. Overall this scheme 
offers rather limited scope for cooling, but as there is consid- 
erable work underway in experiments to reduce the tempera- 
ture of polar molecules to the degenerate regime, this scheme 
might be a useful final stage in the cooling sequence. 



5. Heat capacity 



The heat capacity can be evaluated as 



C 



dE 
df 



T 



N 



dS_ 
df 



(16) 



A' 



i.e. either through the energy functional (9) or the entropy 
(15). We have applied both methods and find they agree. 
These definitions require us to take a numerical derivative, 
which we do by finite difference using two solutions that dif- 
fer in temperature by AT = 10 _4 T. Our results for the heat 
capacity are shown in Fig. 6. 

The behavior seen in the heat capacity is easily understood 
from our entropy results [see Fig. 5(b), and Eq. (16)]. In par- 
ticular the crossover of the curves in Fig. 6(a) (i.e. where the 
difference curves cross the horizontal axis in Fig. 6(b)) oc- 
curs at T w 0.2T F where the peak in entropy difference was 
observed in Fig. 5(b). 



B. Dipolar gas coherence 

Developments in experimental techniques (e.g. see [36- 
40]) have enabled the measurement of correlations in ultra- 
cold gases. Such measurements provide useful many-body 
information (e.g. see [41, 42]). 

The first order coherence of the gas is described by the two 
point correlation function 



G«(x,x')H^(xh/;(x')), 



(17) 



(e.g. see [43]) where ip(x.) is the fermionic quantum field op- 
erator. Introducing center of mass [R = j(x + x')] and rela- 
tive [r = x — x'] coordinates, the correlation function relates 
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FIG. 6. (Color online) Heat capacity versus temperature for a dipolar 
Fermi gas. (a) Heat capacity versus temperature, (b) Change in the 
interacting heat capacity from the ideal gas heat capacity. Aspect 
ratios A = 0.1 (blue/dark grey lines), 1 (green/light grey lines) and 
10 (red/grey lines). Interaction parameter Dt = 0.5 (dotted), Dt — 1 
(solid), D t = 2 (dashed). 



directly to the Wigner function as 

G( 1 )(R,r) = y^3W(R,k)< 



zk-r 



(18) 



While it is possible to experimentally measure the first order 
correlation function at two specific points in a trapped gas 
[40], it is useful to consider a volume-averaged correlation 
function that eliminates the dependence on the center of mass 
coordinate (e.g. see [36-39]). Motivated by this we define the 
volume-averaged first order correlation function 



5 (1) (r) 



JdRG«(R,r) 



JdRjn(R+ir)n(R-ir) 



(19) 



which we have normalized by the density (see Eq. (2.20) of 
[43]) so that g (1) (0) = 1. 



Using Hartree-Fock factorization we can also consider 
higher order correlation functions, such as the second order 
expression 

G< 2 > (x, x') = (ft ( x )^t ( X ')^( X ')^,( X )) , (20) 

= n(x)n(x')-|G (1 >(x,x')| 2 , (21) 

which relates to the system density correlations. The normal- 
ized and volume-averaged second order correlation function 
is given by 



5 (2) W = 1 



JdR|G«(R,r)| : 



|dRn(R+|r))i(R- fr) 



(22) 



Results for the g^(r) and g^ 2 \r) are shown in Fig. 7. The 
low temperature case of g^> (r) considered in Fig. 7(a) ex- 
hibits asymmetry: the length scale for the decay of coherence 
is shorter along the z direction than the x (radial) direction. 
For higher temperatures, this feature is washed out by thermal 
fluctuations and the coherence becomes isotropic. We have 
verified that the low temperature asymmetry in the coherence 
arises from exchange interactions. To do this we set <I>,e(x, k) 
to zero in our calculations (but still retaining the direct inter- 
actions, see Sec. Ill C) and found that the correlation functions 
are isotropic. 

The g' 2 ' (r) correlation function reveals the expected anti- 
bunching behavior for fermions (see Fig. 7(b)), i.e. g^ (0) = 
0. We also observe asymmetry in these second order correla- 
tions at low temperatures similar to that seen for the first order 
coherence function. 

To characterize the first order coherence more thoroughly 
we follow [44] and define the coherence length in the in- 
direction (l co h x ) as 



/drJdR|tf)(R,r)|V 
2/dr/dR|G( 1 )(R,r)| 2 



(23) 



and similarly for the other directions (Note: x is a component 
of the relative coordinate). At high temperatures we have the 
limiting (Boltzmann) behavior l co ] x — ► XdB/V&K (indepen- 
dent of direction), with XdB = h/ '\/2irm,kBT. The results 
for the coherence length are shown in Fig. 8, and reveals the 
dependence of the coherence anisotropy on temperature. 
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FIG. 7. (Color online) Normalized volume averaged correlation 
functions for the dipolar Fermi gas. (a) The first order corre- 
lation function (/^(r) and high temperature limit (black curve) 
e -Tkl /A dB q^ yjjg secon( j or( j er correlation function </ 2 ^(r) and 

high temperature limit 1 — e~ 27r ' r ' AdB . Results given along the x 
axis (dashed) and along the z axis (dotted) for aspect ratio A = 1 
(other aspect ratios show similar behaviour) and Dt — 1. Re- 
sults for T = 0.01T£ (blue/grey) and T = 0.5T£ (red/light grey). 
T = 0.5Tp is used for the high temperature limits. Results for both 
aspect ratios are indistinguishable from their high temperature limit 
at T = T%. 



C. Hartree results 

The results presented so far are computationally demanding 
and it is desirable to produce a simpler theory. We have ob- 
served that, except in the nearly spherical traps, the exchange 
interaction is typically much smaller than the direct interac- 
tion. Thus neglecting the exchange (Fock) term would seem 
to be a reasonable approximation in anisotropic traps, and we 
can use the simplified Hartree dispersion 



e(x,k) 



2 7,2 



Vk 

2m 



[/(x) + $ D (x). 



(24) 



Since the k dependence is of a simple isotropic form it can be 
integrated out to obtain 



n(x) = 



dk 



(27T) 

3-C3/2 ( e 



^dB 



\j*-V*t(x)]/k B T 



(25) 
(26) 



where C,~{z) — Yuk=\(~ -0 ^ /k u is the Fermi function, 
and V e s (x) = U(x.) + $_d(x) is the effective potential. Thus 
the spatial density can be obtained self-consistently using 




Eq. (26), without needing to construct W(x, k) or e(x, k). 



FIG. 8. (Color online) Coherence length of a dipolar Fermi gas as 
a function of temperature: l co h,x (dashed line), l co h,z (dotted line). 
Results are given for D t = 1 and aspect ratio A = 1. Aspect ratios 
A = 0.1 and A = 10 are similar. The black curve shows Xcib/V&k 
for reference. 



To obtain some information about the applicability of the 
Hartree theory, we compare its predictions to the full Hartree- 
Fock calculations for a range of thermodynamic parameters in 
Table II. We observe that agreement between the two theories 
is best for highly anisotropic traps, small interaction strength, 
and high temperatures. Of course we note that in all cases 
the Hartree theory would predict a = 1 (since momentum 
distortion is solely from exchange interactions) and that the 
coherence properties of the gas would be isotropic. 



T = 0.01T 



T = 0.5T£ 



A D t 



^/k B T F 
HF H 



HF H 



E D /Nhuj 
HF H 



fx/k B T F 
HF H 



HF H 



E D /Nfrw 
HF H 



0.1 0.5 
1 



0.945 
0.874 



0.947 
0.887 



0.988 
0.970 



0.989 
0.972 



-3.11 
-7.50 



■3.09 
■7.20 



0.194 
0.168 



0.194 
0.169 



0.998 
0.996 



0.998 
0.996 



-1.07 
-2.26 



-1.07 

-2.25 



0.25 0.5 
1 



0.954 
0.896 



0.957 
0.908 



0.966 
0.921 



0.966 
0.924 



-2.54 
-6.02 



■2.53 
■5.80 



0.198 
0.177 



0.199 
0.178 



0.994 
0.988 



0.994 
0.988 



-0.88 
-1.84 



-0.88 
-1.83 



0.5 0.5 
1 



0.971 
0.931 



0.973 
0.941 



0.941 

0.875 



0.941 
0.880 



-1.61 
-3.86 



■1.60 
■3.72 



0.206 
0.192 



0.206 
0.193 



0.990 
0.979 



0.990 
0.979 



-0.55 
-1.14 



-0.55 
-1.14 



1 0.5 

1 



0.996 
0.983 



0.998 
0.992 



0.922 
0.848 



0.923 
0.853 



-0.20 
-0.82 



■0.19 
■0.77 



0.217 
0.216 



0.218 
0.217 



0.986 
0.972 



0.986 
0.972 



-0.02 
-0.07 



-0.02 
-0.07 



2 0.5 

1 



1.025 
1.043 



1.027 
1.050 



0.923 
0.859 



0.924 
0.863 



1.44 
2.49 



1.44 

2.47 



0.231 
0.243 



0.232 
0.244 



0.985 
0.972 



0.985 
0.972 



0.60 
1.15 



0.60 
1.15 



4 0.5 

1 



1.052 
1.094 



1.053 
1.100 



0.942 
0.898 



0.942 
0.900 



2.85 
5.16 



2.84 
5.10 



0.244 
0.268 



0.244 
0.269 



0.989 
0.979 



0.989 
0.979 



1.15 
2.20 



1.15 

2.20 



10 0.5 

1 



1.074 
1.136 



1.076 
1.141 



0.969 
0.949 



0.970 
0.950 



4.01 
7.20 



4.00 
7.12 



0.255 
0.288 



0.255 
0.289 



0.994 
0.989 



0.994 
0.989 



1.62 
3.08 



1.62 
3.08 



TABLE II. Comparison of Hartree (H) and Hartree-Fock (HF) predictions for thermodynamic parameters. 
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IV. CONCLUSIONS AND OUTLOOK 



In this paper we have developed a finite temperature mean- 
field treatment of a single component trapped Fermi gas with 
dipole-dipole interactions. We have discussed the details of 



our numerical implementation to assist others in implement- 
ing this theory, and have compared to the numerical results 
of Zhang and Yi [29] and variational theories to validate our 
calculations. 

We have considered a range of thermodynamic proper- 



ties of the system. We have found that the long-ranged and 
anisotropic interaction causes many system properties to de- 
pend on the trap geometry, beyond simple scaling with the 
(geometric) mean trap frequency. Interestingly we find that 
dipolar interactions allow the system to be cooled as the har- 
monic confinement is used to squeeze the system towards an 
oblate geometry. 

We have constructed the first and second order correlation 
functions within the Hartree-Fock approximation, and have 
evaluated these in the volume-averaged form relevant to ex- 
periments. We find that the exchange interactions give rise to 
anisotropy in these correlation functions, clearly revealed in 
the coherence length of the system. 

Importantly, our results show the sensitivity of many of the 
features of the dipolar Fermi gas to temperature. In particular, 
over a wide range of parameters and observables we find that 
the strongest effects of the dipolar interactions are apparent for 
T < Q.2Tp. Additionally, we find that many of the qualitative 
predictions made for the uniform gas are strongly modified in 
the trapped system by strong direct interaction effects, partic- 
ularly in oblate traps. Finally, we have investigated a simpler 
Hartree theory that avoids the technical complexities associ- 
ated with evaluating the exchange interactions, and is much 
faster to calculate. We have validated its applicability and ac- 
curacy for studying thermodynamic properties in anisotropic 
trapping potentials where the exchange effects are less signif- 
icant. 

With current fermionic polar molecule experiments pro- 
gressing towards degeneracy, and suggestions that very flat 
traps will be necessary to reduce inelastic collision processes 
[3], these results will be useful to understand the system be- 
havior and to perform accurate thermometry. 
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of their known symmetry. Additionally, the grid allows us to 
implement numerical quadrature, which in the standard form 

is 



N 



dxw(x)f(x) 



3=1 



Wjf{Xj), 



(A3) 



where the roots ({xj}), weight function (w(x)) and weights 
(wj) define the quadrature. In the following subsections we 
introduce the two types of numerical grids we use in detail 
and discuss their related quadrature. 



a. Axial grids 

In the axial direction we use the uniformly spaced position 
grid (i.e. the well-known grid for the discrete Fourier trans- 
form) 



1 



3 



A~ 



3 



1, 



,N, 



(A4) 



where A 2 is the point-spacing. We restrict this grid to positive 
values because of the assumed even symmetry of the functions 
about z = 0. We take the grid for the Fourier transformed 
variable to be 



FzJ ~ tf.A, 



] ~\ 



j = l,...,N z 



(A5) 



Using the trapezium rule with the discrete Fourier grid has 
been shown to be a Gauss-Chebychev quadrature of the first 
kind [45], i.e. quadrature with weight function w(z) = 1 and 
weights w\ = A 2 . A similar quadrature exists in F z space 

withiuf^ = n/N z A z . 

We now consider how to perform the Fourier transforma- 
tion. As the function to be transformed is even we need only 
implement a cosine transform, i.e. 



(A6) 
(A7) 



f(F z )= / dzcos(F z z)f(z), 

J — OO 



Appendix A: Numerical methods 

1. Choice and use of computational grids: underlying 
quadratures 

We represent the Wigner function and the Hartree-Fock dis- 
persion as four dimensional arrays given by 



Wqrst = W(p q ,Z r ,k pg ,k Zt ), 
tqrst — £\Pq i %r i &p s i ^z± ) • 



(Al) 

(A2) 



We have chosen the numerical grids {p q , z r , k Ps , k Zt } to pro- 
vide an efficient representation of these functions making use 



where we have introduced the transformation matrix 

i\ ( r 



T, 



2A 2 cos 



n 



(A8) 



To obtain (A7) we have used the quadrature rule, and assumed 
even symmetry of the function f(z). Similarly the inverse 
transform from f(F z j) — > f(zi) is given by the matrix 



{T- l )> 



N,A, 






i 

l ~2 



1 
] -~2 



(A9) 



Finally, we note that while we have cast this argument in 
terms of performing quadrature, we of course have obtained 
the standard discrete cosine transformation. Thus, to within 



10 



constant prefactors, the forward and inverse transform are the 
so called DCT-IV variant of the discrete cosine transformation 
available in some standard FFT libraries (e.g. see [46]). 

For the case of the axial momentum variable k z , the grid is 
similarly given by 



3-\) 



l,...,JVfc. , 



(A10) 



where A kz is the spacing between momentum points. The 
above discussion regarding the quadrature of quantities on the 
position grid immediately applies to the momentum case with 
the replacements N z — > N&, and A z — > Afc z . 



b. Radial grids 

In the radial direction it is convenient to introduce grids 
based on the roots of the Bessel function to allow an effi- 
cient Hankel transform (i.e. radial Fourier transformation) to 
be performed (also see [34]). For the position variable we use 
a numerical grid of N p points covering the region (0, R) given 
by 



Pj 



ajR 



j = l,...,N p , 



(All) 



where ctj is the jth root of the Bessel function J${z). We take 
the respective grid for the Fourier transformed variables to be 



F - = 21 

p - J R ' 



j = l,...,N p 



(A12) 



Associated with these grids are quadratures with weight 
functions w(p) = p and w(F p ) = F p in direct and Fourier 
space, respectively, with the corresponding weights 



w 



,(i>) 



(J> 



2R 2 



1 



a N +1 Ji{a f 



R 2 Ji(aj) 2 ' 



(A13) 
(A14) 



see [47, 48]. 

The in-plane (i.e. two dimensional) Fourier transformation 
of a radially symmetric function, f(p), is given by 



f{F p ) 



10 

2tt 



'"dcjie^^fip), 



dppJ (F p p)f(p) 



(A15) 
(A16) 



where 4> is the in-plane angular coordinate. The term in the 
square brackets in Eq. (A16) is a Hankel transformation. Mak- 
ing use of the quadrature (A16) is approximately given by 



3 

where the transformation matrix is 

4ttR 2 Jo (a-iOj/a-iVp+i) 



H,. 



a 



N p 



Ji{* 3 ) 2 



(A17) 



(A18) 



The inverse transformation, i.e. f(F p ,i) —> f(pj), is given by 
the matrix 



H- 



1 J {a-iUj/uN p +i) 



'v AttR 2 



Ji(a 3 ) 2 



(A19) 



For the case of the radial momentum variable k p , the grid 
is similarly given by 



1 Pj 



a N kn +l 



J = 1, 



,N k 



(A20) 



where K is the range of the momentum grid. The above dis- 
cussion regarding the quadrature of quantities on the position 
grid immediately applies to the momentum case with the re- 
placement N p — s- Nk and R —J- K. For future reference, we 

denote w\ as the quadrature weights. 



2. Evaluating the direct term 

Here we demonstrate how (5) is evaluated numerically. 
First, we obtain the density, 

n(p,z)= dk z dk p ^W(p,z,k p ,k z ), (A21) 



/o ./o 



2ir 2 



(A22) 



si 



where we have made use of the momentum variable quadra- 
tures to perform the integrations. 

Then we note that the convolution is given by Fourier trans- 
forms as $d(x) = J r_1 [?7dd(k)n(k)]. Ronenetal. [34] have 
shown that the discontinuity of t/dd(k) at the origin results in 
slow convergence without a cutoff in the real space interaction 
potential. Using a cutoff, L, that is larger than the system has 
no physical consequences, and results in the corrected Fourier 
transformed interaction: 



^(k) 



1 



:C, 



<l,l 



cos(Lk) sin(Lfc) 



{Lkf 



{Lkf 



(3cos0£-l), 
(A23) 



Then we have that (5) is given by 

®D,mn=®D{Pm,Zn), (A24) 

= T.i H ^) n A T ^) n] ^ d {F p „F z ,) 



1.1 
X 



7 fliq-L jr'l'qr ■ 



(A25) 



qr 



Our calculations for anisotropic traps require the cutoff to 
be larger than any of the dimensions of our cloud. As our 
aspect ratios are up to a factor of ten, it is computationally 
efficient to avoid evaluation of the terms in (A25) where the 
density is known to be negligibly small. 
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3. Evaluating the exchange term 



We now consider (6). Unlike the direct term, which acts on 
operate on a marginal (i.e. the density), the exchange interac- 
tion acts on the complete Wigner distribution. We adopt the 
same notation for the direction term, i.e. 



<I> 



E,qrmn 



^E{Pq,Z r ,k Pri 



(A26) 



where we have integrated out k', and used that W(x, k) is 
even in k z and has no kj, dependence, and defined 



A(fcp, fc z , ftp, k z ) 



Cdd 

3(2tt) 2 X 

Z(k z -k' z ) 2 



(k 2 + k' 2 - 2k z k' z ) 2 - Ak 2 p k 12 



(A30) 



We note that Eq. (A29) is consistent with the functional 
derivative of the expression Zhang and Yi have developed for 
the exchange energy (see Sec. IV of [18]). 

Numerically the exchange term is evaluated as 



$ 



E,qrmn 



A^^W^WqrstAmns 



(A31) 



where A mnst = A(k Pm , k Zn ,k Ps , k Zt ). 



We consider an analytic simplification of the exchange term 



4. Choice of number of grid points 



$ B (x,k) 



dk' 

(2^)3 



Ty(x,k')f/dd(k-k'), (A27) 

3(fc z - k' z f 



C dd f dk' 

3 J (2tt) 3 v 



|k-k' 



dk' 



k' p dk p W(x,k p ,k' z )x 



1 , 
(A28) 



[A(k p , k z ,k' p , k' z ) + A(k p , k Zl k' p , -k' z )] , (A29) 



Using the gaussian ansatz of [28], we compared the total di- 
rect and exchange energies to the exact result at T — 0.5TJL 
We found that the relative error in the direct energy decreased 
exponentially with {N p ,N z } and was less than 10 -12 with 
{N p ,N z } = {24,24}. The relative error in the exchange 
energy decreased algebraically with {N k , Nk, }, Nk = 
0.6N kz was optimum and with {N kp ,N k J = {48,80} the 
relative error was up to 1CP 4 . 

For the results presented in Section III, the Wigner func- 
tion and dispersion relation are calculated on four dimensional 
grids with sizes of {N p ,N z ,N kfl ,N k J = {40,40,48,80}. 
For each temperature these calculations can take approxi- 
mately 10 hours to converge to self-consistency on an eight- 
core 2.8 GHz Xeon processor. 
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